root.dir <- here::here()
knitr::opts_chunk$set(echo = TRUE, root.dir=root.dir)
knitr::opts_knit$set(root.dir = root.dir)
knitr::opts_chunk$set(message = FALSE)
knitr::opts_chunk$set(warning = FALSE)
The data/haplotype_amino.fas from haplotype analysis is used for generating haplotype tree and assessing correlation between haplotypes.
library(adegenet)
library(ips)
library(ggplot2)
library(ggtree)
library(heatmaply)
nbin<-fasta2DNAbin("data/haplotypes_AA_names.fas")
an<-as.alignment(nbin)
nm<-as.matrix(an)
nbinmat<-as.matrix(labels(nbin))
class(nbin)
dnbin<-dist.dna(nbin, model = "K80")
tree<-njs(dnbin)
ggt<-ggtree(tree, cex = 0.8, aes(color=branch.length))+scale_color_continuous(high='lightskyblue1',low='coral4')+geom_tiplab(align=TRUE, size=5)+geom_treescale(y = -1, color = "coral4", fontsize = 7)
njmsaplot<-msaplot(ggt, nbin, offset = 0.009, width=1, height = 0.5, color = c(rep("rosybrown", 1), rep("sienna1", 1), rep("lightgoldenrod1", 1), rep("lightskyblue1", 1)))
njmsaplot
pdf("figures/haplotype_tree.pdf", width = 11, height = 9, paper = "a4", pointsize=15)#save as pdf file
njmsaplot
dev.off()
sat2 <- NULL
for (i in 1:nrow(nm)) {
sat2[i] <- paste(nm[i, ], collapse="")
}
sat2 <- toupper(sat2)
sat3 <- unique(sat2)
comat = matrix(nrow=length(sat3), ncol=length(sat3))
for (i in 1:length(sat3)) {
si <- sat3[i]
for (j in 1:length(sat3)) {
sj <- sat3[j]
difcnt = 0
s1 = as.vector(strsplit(as.character(si), ""))
s2 = as.vector(strsplit(as.character(sj), ""))
for (k in 1:length(s1[[1]])) {
if (s1[[1]][k] != s2[[1]][k]) {
difcnt = difcnt + 1
}
comat[i, j] = difcnt
#print(paste(i, " ", j, " ", difcnt))
}
}
}
#comat is Hamming distance matrix
colnames(comat)<-nbinmat
heatmap<-heatmaply_cor(cor(comat), file= "figures/heatmap.pdf", xlab= "haplotypes", ylab="haplotypes",
k_col=3,
k_row=3, margins =c(10,10,15), fontsize_row = 8,
fontsize_col = 8)
heatmap
Figures are available haplotypes_trees and Haplotype_correlations